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Abstract 

Verification of temporal logic properties plays a crucial role in proving the desired behaviors of hybrid sys¬ 
tems. In this paper, we propose an interval method for verifying the properties described by a bounded 
linear temporal logic. We relax the problem to allow outputting an inconclusive result when verification 
process cannot succeed with a prescribed precision, and present an efficient and rigorous monitoring algo¬ 
rithm that demonstrates that the problem is decidable. This algorithm performs a forward simulation of a 
hybrid automaton, detects a set of time intervals in which the atomic propositions hold, and validates the 
property by propagating the time intervals. A continuous state at a certain time computed in each step is 
enclosed by an interval vector that is proven to contain a unique solution. In the experiments, we show that 
the proposed method provides a useful tool for formal analysis of nonlinear and complex hybrid systems. 

Keywords: Hybrid systems, interval analysis, linear temporal logic, bounded model checking. 


1 Introduction 

Reasoning of the temporal logic properties in a hybrid system is a challenging and 
important task that lies in the intersection among computer science, numerical anal¬ 
ysis, and control theory. Various methods for falsification of hybrid systems with 
temporal properties have been developed, e.g., [21,20,6,27], and these methods en¬ 
able verification of various properties (e.g., safety, stability, and robustness) of large 
and complex systems. The state-of-the-art tools are based on numerical simula¬ 
tions whose numerical errors often produce a qualitatively wrong result and become 
problematic even in a statistical evaluation. 

A fundamental process in formal methods for hybrid systems is computation of 
rigorously approximated reachable sets. The techniques based on interval analysis 
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(Section 3) have shown practicality in the reachability analysis of nonlinear and 
complex hybrid systems [8,5,24,15,3,11,12], In these frameworks, computation is 5- 
complete [10]: function values are allowed to be perturbed within a predefined 6 G 
M>o, and, by setting bounds in a problem description, many generically undecidable 
problems become decidable. The d-complete verification of generic properties other 
than reachability is a challenging topic. 

In this paper, we present an interval method for verifying bounded linear tem¬ 
poral logic (BLTL) properties (Section 5) for a class of hybrid automata (Section 4). 
Our method computes three values in a reliable manner: the algorithm assures the 
soundness using interval analysis when the result valid or unsat is output; other¬ 
wise, unknown is output when the verification process reaches a prescribed precision 
threshold. We present an algorithm (Section 6) based on the forward simulation of 
a system. It encloses a trajectory with a set of boxes (i.e., interval vectors) and also 
ensures the unique existence property (i.e., we ensure that a unique state is enclosed 
in a box corresponding to each initial value) for each step of the simulation. For 
each atomic proposition involved in a property ip to verify, the algorithm obtains 
an inner and outer approximation of the time intervals in which the proposition 
holds. Next, the set of time intervals is modified according to the syntax of the 
property ip, and finally the algorithm checks whether ip holds at the initial time. 
Using our implementation, we show that nonlinear models are verified and the nu¬ 
merical robustness of a trajectory is assured (Section 7). Although our method is 
simple, it enables reliable analysis of a set of trajectories and provides a foundation 
for validated model checking and controller synthesis. 


2 Related Work 

Many previous studies have applied interval methods to reachability analysis of 
hybrid systems [8,5,24,15,3,11,12]. The outcome of these methods is an over¬ 
approximation of a set of reachable states with a set of boxes. In interval analysis, 
a computation often provides a proof of unique existence of a solution within a 
resulting interval. This technique also applies in interval-based reachability analy¬ 
sis [15,14], but it is not considered in most of the methods for hybrid systems. Our 
method enforces the use of the proof to verify more generic temporal properties. 

Reasoning of real-time temporal logic has been a research topic of interest [2,25]. 
Numerical method for falsification of a temporal property is straightforward [16]. 
It simulates a trajectory of a bounded length and checks the satisfiability of the 
negation of the property described by a bounded temporal logic. This paper presents 
an interval extension of this falsification method. 

A tree-search method for searching witness trajectories [21], a falsification 
method based on a Monte-Carlo optimization technique [20], and statistical model 
checking methods [6,27] have been proposed. These methods have been shown their 
practicality in the verification of realistic nonlinear models; however, their imple¬ 
mentations are based on numerical simulations and might suffer from numerical 
error. Applications of our interval method include an integration with these statis¬ 
tical methods to achieve both reliability and practicality. An integrated statistical 
and interval method was also proposed in [26] for reachability analysis. 
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Notions of robustness have been proposed to facilitate the simnlation-based ver¬ 
ification of temporal properties [9,7,20]. In these works, the degree of robustness 
is represented as a distance between a trajectory and a region where a proposition 
holds. A non-robust trajectory, which is computed numerically, is likely to be in¬ 
consistent with the considered model due to numerical errors. Our method ensures 
a robustness rigorously by verifying that a trajectory intersects with each boundary 
in the state space. 

There exist a few methods for model checking of temporal logic properties [23,4]. 
[23] proposed a method specialized in stability properties, which is described as a 
specific form of temporal logic formula. [4] proposed a method that translates a 
verification problem into a reachability problem with the /c-Liveness scheme, which 
is incomplete in general settings. Our method can be viewed as a bounded model 
checking method that validates a bounded temporal property by ensuring that all 
trajectories that emerge from an initial interval value satisfy the property. 


3 Interval Analysis 

This section introduces selected topics and techniques based on interval analy¬ 
sis [17,19]. The techniques are used in the proposed method in Section 6. 


3.1 Basic Notions and Techniques 

A (bounded) interval a = [a, a] is a connected set of real numbers {6 G M j a < 
6 < a} and I denotes the set of intervals. For an interval a, a and a denote the 
lower and upper bounds; the width is defined as a —a] and int a denotes the interior 
{5 G M j a < 6 < a}, [a] denotes a point interval [a, a]. For intervals a and b, 
d{a,b) denotes the hypermetric between the two, i.e., max(ja — b\,\a — b\). For 
a set S' C M, DS denotes the interval [inf S, sup S]. All of these definitions are 
naturally extended to interval vectors; an n-dimensional box (or interval vector) 
a is a tuple of n intervals (ai,...,a„), and denotes the set of n-dimensional 
boxes. For a G and a G I”", we use the notation a £ a, which is interpreted as 
VzG{ 1, ..., n} Oi G ttj. In an actual implementation, the bounds of intervals should 
be machine-representable floating-point numbers and other real values are rounded 
in the appropriate directions. 

For a function / : M"' —)■ M, f : —)• I is known as an interval extension of / 
if and only if it satisfies the containment condition Va G Va G a (/(a) G /(a))- 
This definition is generalized to function vectors f : M” —)> where Uf G N>i. 
Given intervals a, 6 G I, interval extensions of four operators o G can 

be computed as nla o b, aob,a o b,a ob} (we assume 0 0 6 for division). 

For arbitrary intervals a,b,d G I, the extended division Did £ d \ 3a G a 36 G 
b a = bd} can be implemented as follows (see Section 4.3 of [19]): 


ExtDiv(a, b, d) 


a/b n d 

□ (d\ (a/6, a/6)) 

□ (d\ {d/b,d/b)) 

d 

\ 


if 0 0 6 
if a > 0 G 6 
if a < 0 G 6 
if 0 G a, 6 
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In the second and third cases, when 6 = 0 (resp. b = 0), we set a/6 and a/6 as —oo 
and oo (resp. a/6 and a/6 as oo and —oo). 

Given a differentiable function /(a) : M —?• M and a domain interval a, a root 
a G a of / such that /(a) = 0 is included in the result of an interval Newton operator 

an (a - f{a)/£f{a)) , 

where a ^ a, and f and are interval extensions of / and the derivative of /. 
Iterative applications of the operator will converge. Let a' be the result of applying 
the operator to a. If a' C int a holds, a unique root exists in a'. 

3.2 ODE Integration 

An initial value problem (IVP) for an ordinary differential equation (ODE) is spec- 
ihed by a triple (to, .^) consisting of an initial value xq G at time to gR. and 
a flow function F : M” —)• M” (assume Lipschitz continuity). Given a time interval 
t G I and a continuous trajectory x{t) : t —)■ M”, the satisfaction for IVP-ODEs is 
defined as 


x,t \= {to,xo,F) iff x{to) = xo/\yiGt -^x{i) = F{x{i)). 

Given to ^ I ®o G I", we can consider a parametric IVP-ODE (to, xq, F), where 
the initial condition is parameterized, and its satisfaction relation is defined as 

x,t \= {to,xo,F) iff 3toGto BxqGxo x,t\= {to,xo,F). 

TSt{tQ,XQ, F) denotes the set of satisfied trajectories on t. 

Using the tools based on the interval Taylor methods, e.g., GAPD'^ and VN¬ 
ODE [18], we can obtain an interval extension X : I —)• I"" of solution trajectories in 
TSt{to,xo, F). Given t' G I, such tools compute a value X(t') by performing the 
stepwise integration of the flow function F from the initial time to to time t'. In 
the stepwise computation of the interval Taylor methods, the unique existence of a 
solution is verified for a box enclosure computed in each step based on the Picard- 
Lindeldf operator and Banach’s fixpoint theorem. Accordingly, when an interval 
enclosure X(t') (assume il > to) is computed with an interval Taylor method, the 
following property holds: 

Vto Gto Vxo G®o dunique x G (t^ —?■ X(T)) x, t' |= (to, a^o, F)- 

In principle, if F is Lipschitz continuous and we can assume an arbitrary pre¬ 
cision, we obtain an arbitrary narrow interval enclosure X([t]) for t G M. However, 
since the implementations use machine-representable real numbers, it may fail to 
compute an enclosure in the process that verifies the unique existence property, 
even with the smallest step size. 
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4 Hybrid Automata 

We model a hybrid system as a hybrid automaton [1], For simplicity in this paper, 
we consider deterministic systems, i.e., the location invariant is the negation of 
guard conditions and two guards do not overlap in a location. The proposed method 
can be extended to handle non-deterministic systems, e.g., by enumerating possible 
paths and computing a trajectory enclosure for each path. 

Definition 4.1 A hybrid automaton is a septet 

HA := (^Q, X, X, Init, {Fq}q£Q, {Gq^q'}q£Q^q'£Q, {Rq^qi}q£Q^qi£Q'j , 

that consists of the following components: 

• A hnite set of locations Q = {qi,..., q-nq}- 

• A vector of real-valued variables x = (xi,..., Xn)- 

• A domain A C I” for the valuation of the variables. 

• A set of initial values Init C {q}xX where q € Q. 

• A set of vector fields Fq : X ^ X (assume Lipschitz continuity). 

• A set of guards Gq^qi C X described by a condition of the form g{x) = 0 A h{x) < 0 
where g,h : X ^ M.. 

• A set of reset functions Rq^qi : A —)> A. 

Behaviors of the states a G QxX over the timeline are formalized as trajectories. 
In this work, we assume that there are no consecutive multiple discrete changes. 

Definition 4.2 Given an HA, an initial state {qo,SQ) G Init, and a time interval 
t = [0,tmax] (imax £ lK>o)j & State at each time t G f is determined as a trajectory 
(q, s), which consists of a location trajectory q : t —)• Q a continuous state 
trajectory s : t —)• A. The value of the trajectory is dehned recursively as follows: 

(q(0),s(0)) := iqo,So), 

(q(t),s(t)) := a G QxX s.t. 3t' G{0,t) 3a' gQxX (q(t'), s(t')) a' a, 

where the relation cti -A cj 2 G (QxX) xtx (QxX) is given by the following rules: 

/(O) = s Vt G [0, t] Af{i) = Fq{f{i)) 

Gq,qfis) Rg^qfis) = s' _ VtG[0,t) Vg'G Q ^Gq^qfifji)) 

(q,s) ^ (q',s') iq,s){q,f{t)) 

Note that the second rule also applies for t = 0. The set of trajectories of length 
imax is denoted by TSt^,^^{HA). 

When a discrete change is applied at time t, the state (q(t),s(t)) overwrites 
the state a' before the discrete change. An HA has a unique trajectory (q, s) starting 
from an initial state {qo,so) G Init because we have assumed that two guards do 
not hold simultaneously. 
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frtve U[o,5]Xi>2 i i i i 

-.(frue U[o,5]Xi>2) - 

frtve U[o,io]-'(frue U[o,5]Xi>2) -■ = 

^{true U[o,io] ^{true U[o,5]>fi>2)) i - 

Figure 1. Verification of the bouncing ball example 

Example 4.3 We model a bouncing ball on a moving table as an HA\ 


X := (xi, X2, X 3 ) e X := ([-1,10], [-10,10], [0,1000]) 

L := {g} 

Init :={g}x([2,7],[0],[0]) 

Fq := (x 2 , — 1 + 0.04x2 sgnx 2 ,1) 

Gq^q := xi — sinxs = 0 A X 2 — COSX 3 < 0 
Rq,q := (xi, — 0 . 9 X 2 + 1.9 cos X 3 , X 3 ) 

Variables xi, X 2 , and X 3 represent the height and velocity of the ball, and the (global) 
time, respectively. Air resistance is considered in the dynamics. The height of the 
table sinusoidally oscillates within [—1,1] and is represented as sinx 3 . The second 
proposition of the guard is to forbid the guard to hold right after a discrete change. 
Possible trajectories of xi and X 2 are illustrated in Figure 1. 


5 Bounded Linear Temporal Logic 

We consider a fragment [16] of the real-time metric temporal logic [2] such that the 
temporal modalities are bounded by an interval t = [t, t] such that the bounds t, t 
are in Q. We refer to the logic bounded linear temporal logic (BLTL) as in [27]. 

Definition 5.1 We consider constraints in the real domain as atomic propositions. 
The syntax of the BLTL formulae is defined by the grammar 

(/? ::= true j p j V </? j -■(/? j ip Ut ip p ::= /(x) < 0 1 /(x) < 0 

where p belongs to a set of atomic propositions AP^,, U* is the “until” operator 
bounded with a non-empty positive time interval t £ I, x = (xi,..., x„) is a vector 
of variables, and / : M” —)• M. We use the standard abbreviations, e.g., ipi A ip 2 := 
V ^ip 2 ), ■= true Ut (“eventually”), and Gtip '■= ^ft^ip (“always”). An 

equation /(x) = 0 can be encoded as /(x) < 0 A — /(x) < 0. 
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5.1 Semantics 

The necessary length ||(/?|| of trajectories for checking a formula ip is inductively 
defined by the structure of the formula: 

IIpII :=0 \\<piy <P 2 \\ := max(||v7i||, ||v? 2 ||) 

lh¥?|| := ||VJ|| W^Pi Utv?2|| := max(||(/?i||, \\(p2\\)+i 

A map O : AP^ —)■ 2^ corresponds each proposition p G to a set 0{p) = {sG 
X I p(s)}. Let (q,s) be a trajectory in TS\\^p\\{PlA) and (/? be a BLTL property. We 
have a satisfaction relation defined as follows: 


s, t = true 


s, t 1= p 

iff s{t) G 0{p) 

S, t 1= VPl V ip2 

iff s, t \= ipi V s,t = ip 2 

S, t 1= ^if 

iff s, t ^ (/p 

S, t 1= ifi Ut ip2 

iff 3t' G (t -|- t) s, t' = ip 2 A fit" G [t, t'] s, t" = ipi) 

ip 2 intuitively means 

that (assuming we are at time t) ip 2 will hold within the 


time interval t+t and ipi always hold until then. We also have a validation relation 
defined as: 

HA^p> iff V(q, s) G T5||^||(//A) s, 0 h VP 
5.2 Method for Monitoring BLTL Formulae 

Our interval method is based on the method proposed in [16] that decides whether a 
trajectory satisfies a BLTL property. In this section, we explain this basic method. 
First, we introduce the notion of consistent time intervals against BLTL formulae. 

Definition 5.2 Let (q, s) be a trajectory of length tmax and be a BLTL formula. 
We say that a left-closed and right-open interval [t,t) C M>o is consistent with ip 

iff VtG [t,t) s,t \= T- 

The satisfiability of a property (phy a. trajectory is checked as follows: 

(i) For each atomic proposition p in (/?, monitor the trajectory of length ||(v>|| and 
identify a non-overlapping set of consistent time intervals Tp = {fi,..., fnp}- 

(ii) Following the parse tree of (/p in a bottom-up fashion, compute a set of consistent 
time intervals of ip. For each construct of BLTL, compute as follows: 

r-,<p := M>o \ Tp T^\yiP2 ^ T^2 

'^v’iUtV 2 ■= {Shiftt(fi n t 2 ) n I ti G r^i,i 2 e ^^ 2 } 

where Shiftt(s) := [s — t,s — t)ri M>o. 

(hi) Check whether the smallest time interval in T^p contains time 0. If yes, ip is 
satisfied; otherwise, it is not satisfied. 

Example 5.3 We verify the property 

G[o,io]F[o, 5] 2 - xi<0 = ^(trueU[o,io]^(trueU[o,5]2-xi<0)) 
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for the model in Example 4.3. Computation with the monitoring method (which is 
extended to an interval method) is illustrated in Figure 1. 

6 Interval-Based Simulation and Monitoring Method 

In this section, we propose an interval extension of the monitoring method in Sec¬ 
tion 5.2. 

In Step (i) of the method, given an HA and a BLTL property (p, we first simulate 
the HA for length ||(^||. Our simulation method computes an over-approximation of 
a trajectory of HA in which the existence of a unique trajectory is verified, i.e., it 
verifies the property 

\/{qo,xo)£lnit dunique (q,s) € TS'||,^||(i^^) q(0) = go A s(0) = xq 

meaning that for each initial value, there exists a unique trajectory of length ||(^||. 

Next, our method identifies the time intervals that are consistent with an atomic 
proposition in Lp (Definition 5.2). A consistent time interval [t,t) is, in general, not 
representable in an actual implementation; therefore, we approximate it by a pair 
of intervals u and u' such that each of them encloses the boundaries t or t. Given 
an atomic proposition f{x) o 0 (o E {<, <}) and a trajectory s, we search for a 
boundary enclosure u such that f{s{u)) 3 0 where / is an interval extension of /. 

In Step (ii), the set of consistent time intervals is updated to be consistent 
with (p. This computation requires u to contain a unique boundary point, and 
thus a naive over-approximation u does not suffice because an interval extension 
f{u) may contain 0 although f{u) does not contain a boundary or contains several 
boundaries. Thanks to interval techniques, our method verifies the unique existence 
of a boundary in u. Finally, as an extension of the above property, our method 
verifies 

'^{qo,xo)£ Init dunique (q, s) G TSj^j{HA) q(0) = go A s(0) = xq A s, 0 ^ 

The proposed method has some limitations. First, it is a semi-decision pro¬ 
cedure that may output an inconclusive result (unknown) because of a failure in 
the verification of unique existence; both the procedures for enclosing a continuous 
trajectory and enclosing a time where a discrete change occurs may cause errors. 
However, this mechanism is valuable in terms of reliability and complexity of the 
problem; a non-robust trajectory and a zeno HA will be rejected as an error in 
the verification process. In practice, when addressing a nonlinear HAs, the method 
may only work successfully with a sufficiently small subset Init' C Init of initial 
values. In this way, the method can be still used for sat/invalid checking. Second, 
the method is a bounded model-checking method in the sense that the domain X of 
the variables is bounded, and it assumes a bounded length and a number of discrete 
changes in a trajectory. 
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Input: HA, 

Output: valid, unsat, or unknown 

1 : t ;= 0; q:=qo; x := xq; T:={0, ...,0} 

2 : while t < ||(/?|| do try 
3: tc := M 

4: for q' ^ Q do {Find zeros for each edge} 

5: t'c := SearchZero(X,Fg,Gg,q/, [t, \\ip\\]) 

6: if t'c / 0 A < tc then 

7: q” ■.= q'] tc := 

8 : else if > tc then 

9: return unknown 

10 : end if 

11 : end for 

12 : for p = f o 0 £ AP^p do {Find boundaries of APs} 

13: t'^-.= [t,tc] 

14: loop 

15: tc •= SearchZero(X, Fq, / = 0, t'^) 

16 : if t'c = 0 then break end if 

17: tc := [tcjtc]; Tp := Tp U {tc} 

18 : end loop 

19: end for 

20 : a; := Jump(X, tc, AP(^, T) {Discrete change} 

21 : q := q”; t := F 

22 : catch error then return unknown end try 
23: end while 

24: return Anallntervals(T) 

Figure 2. Monitor algorithm 


6.1 Main Algorithm 

Given an HA and a BLTL property ip, the proposed Monitor algorithm (Figure 2) 
outputs the following results: valid that implies HA |= (p] unsat that implies HA ^ 
-■(/?; or unknown when the computation is inconclusive. 

An iteration of the outmost loop corresponds to a continuous phase of the tra¬ 
jectory and a discrete change. At Lines 4-11, each guard for a possible transition 
is evaluated. The SearchZero algorithm is described in Section 6.4 and will return 
a time interval within which the guard is satisfied or 0 if no state satisfies the 
guard. Next, the algorithm attempts to decide the earliest time interval by check¬ 
ing whether is strongly less than tc- If two guard crossings are too close, so two 
crossing time intervals overlap, the algorithm returns unknown. At Lines 12-19, for 
each atomic proposition of ip of the form / oO (o £ {<,<}), the algorithm searches 
for a boundary where the sign of / changes. Because several boundaries can exist in 
a continuous phase, the inner loop searches for all of them. The detected time inter¬ 
vals are saved in the set T associated with the atomic propositions. At Lines 20-21, 
the discrete change between the locations q and q” is computed by evaluating an 
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interval extension of i?g^g"(X(tc)). A jump of state might switch the state of an 
atomic proposition; if such a switch exists, the Jump procedure should detect and 
record it in T. Finally, boundary points (which are enclosed by intervals) of the 
consistent time intervals saved in T are analyzed by the Analintervals procedure 
(Line 24, Section 6.2). The procedures X (Section 3.2) and SearchZero (Section 6.4) 
may results in errors. These errors are caught by the catch clause at Line 22. 

6.2 Evaluation of BLTL Properties 

The BLTL evaluation explained in Section 5.2 can be implemented as a rigorously 
approximated procedure Analintervals. 

First, we approximate a set of consistent time intervals = {ti,... by 

= {ui,u']^,... such that G I, G w*, L G Ui < u', and 

E'i < for i G {!,... ,n^p}. T = 0 and T = M>o are approximated as T = 0 
and T = {[0]}, respectively. For Ui, u[ in and a continuous trajectory s, 
Vt G [ui,«') s,t \= if holds. 

Next, the evaluation on the set of time intervals in Step (ii) is extended to 
address the approximated sets. The set of the inverted time intervals for -k/j can 
be represented as T-,,^ = {[0], ui,..., if Ui > 0, and ..., if 

ui = [0]; otherwise, the evaluation results in unknown if ui 3 0. The union of two 
sets of time intervals for ipi V ip 2 can be implemented as a merge and sort process of 
the two approximated sets. The Shift* procedure for ipiyJtP 2 can be implemented 
as translations of the time intervals Ui and ul for t and t, respectively. Some more 
case analyses should be applied, e.g., when the intervals in become redundant 
or when they overlap. 

Finally, we obtain and conclude that ip is valid if ui < 0 < vfp, it is unsat if 
Tip = 0 or 0 < or the satisfaction is unknown if 0 G [ui,ui). 

6.3 Computation of a Continuous Trajectory 

If the system is in a location q € Q and the value of the state variables is a; G I” at 
time t G I, then the subsequent continuous evolution is specified by an IVP-ODE 
{t,x,Fq). Next, we can obtain an interval extension X : I —?■ I”' of the continuous 
trajectories in TS'j* ||(^||](t, F)j) as described in Section 3.2. 

6.4 Evaluation of Boundary Conditions 

Guards and atomic propositions can be treated as boundary conditions in the 
state space X C ME of the form 

B{x) := g{x) = 0 A h{x) < 0, 

where g : X ^ M and /i : A —)• M. We propose the SearchZero algorithm shown 
in Figure 3 for searching the intersection between a trajectory and a boundary 
condition. Inputs to the algorithm consist of an interval extension of the continuous 
trajectory X : I —)• I"", a vector field of the current location F : A —?■ A, the boundary 
condition B(x), and a time interval tjnit G I to be searched. SearchZero searches for 
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Input: X:I —F : X ^ X, g = 0 A h<0, tjnit £ I 
Output: t G I U {0} 

Parameter: e G Q>o, ^ £ (0,1) 

1: t := tinit 

2 : repeat {Lower bound reduction} 

3: ^old •— t 

4: dg := Dt(c/, X, F, t); dh := Dt(/i, X, F, t) 

5: t := t + ExtDiv(—fif(X(t)), dg, t — t) 

6: t := t + ExtDiv(—/i(X(t)) — [0, oo], dh,t — t) 

7: until d{toid,t) < e 

8: if t = 0 then return 0 end if 

9: t := t 

10 : loop {Unique solution existence verification} 

11 : dg := Dt{g,X,F,t) 

12: if dg 3 0 then error end if 

13: t':= t-g{X{t))/dg 

14: if t' C intt then t := t'; break end if 

15: tbak := t 

16 : t := tinit n lnflate(t', 1+0) 

17: if d{t,t') > (1 —0) tbak) then error end if 

18 : end loop 

19: if sup/i(X(t)) > 0 then error end if 
20 : return t 

Figure 3. SearchZero algorithm 

the earliest time interval t C such that the state X(t) encloses a unique solution 
of the boundary condition, i.e., 

t = □{minjtGtinit | B{s{t))} \ se TSt;^^,{to,xo, Fg)] , (1) 

where {tQ,XQ, Fq) denotes the IVP-ODE of the current location. Moreover, 
SearchZero verifies the following property: 

yseTSti^,^{to,xo,Fq) dunique tet B{s{t)) (2) 

Otherwise, SearchZero returns 0 when the boundary condition is unsatisfiable, i.e., 

VSG TSti,,;^{to,Xo,Fq) VfGtinit ^B{s{t)). (3) 

Lemma 6.1 If SearchZero returns a non-empty interval t, the properties (1) and 
(2) hold. If it returns 0, the property (3) holds. 

To justify the soundness, we describe some details of the algorithm. At Lines 2- 
7, the time interval t is filtered repeatedly using an interval Newton operator. At 
Line 4 (and at Line 11), given a function g, the Dt procedure computes an interval 
enclosure of the derivative ^g{s{t)) over the time interval t using the chain rule 

ig{s{t)) = £g{s{t)) . is{t) C £giX{t)).F{X{t)). 
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Next, at Lines 5 and 6, the interval Newton is applied. The extended division 
(Section 3) is used to implement the interval Newton to handle the numerator 
intervals dg,dfi containing zero. Because we expand the interval Newton on the 
lower bound t and the extended division encloses the values in the domain t — t, the 
resulting t is filtered its inconsistent portion, without losing the solutions or being 
expanded. When the interval Newton results in 0, SearchZero also returns 0 to signal 
the unsatisfiability. At Line 9, because t may contain several solutions, t is reset to 
the lower bound as a starting value to compute an enclosure of the earliest solution. 
At Lines 10-18, SearchZero applies the interval Newton method with the inclusion 
test to prove the unique existence of a solution within the contracted interval t'. 
This interval Newton verification is repeated with an inflation process of the time 
interval (see [13] for a detailed implementation). When reaching Line 19 with no 
error, the time interval t is a sharp enclosure of the first zero of g{s{t)) = 0. It 
remains to check that the inequality constraint h{s(t)) < 0 is satisfied inside t. 

When SearchZero is implemented with machine representable real numbers, or 
when there is a tangency between the trajectory and the guard constraint, a com¬ 
putation may result in an error. Line 12 of SearchZero may give rise to error if 
the derivative on an (inflated) time interval contains zero. At Line 17, we limit the 
number of iterations according to whether the inflation ratio reaches the threshold 
as proposed in [13] . 

7 Experiments 

We have implemented the proposed method and experimented on two examples to 
confirm the effectiveness of the method. Experiments were run using a 2.4GHz Intel 
Core i5 processor with 16GB of RAM. 

7.1 Implementation 

We have implemented the proposed algorithms in Figures 2 and 3 in OCaml and 
C/C-I-+. The CARD library was used for solving ODEs. Parameters tmim e, and 9 
should be configured. Each parameter corresponds to the smallest integration step 
size that CAPD can take, the threshold used in Figure 3, or a threshold used in 
Inflate. In the experiments, these parameters were set as fmin := 10“^^, e := 10“^^, 
and 6 := 0.01. 

For most of models, our implementation only accepts small interval values, as 
reported in the next section, because an interval enclosure of the state after a 
continuous evolution and a jump expands by quite a large amount, and thus, the 
verification process in SearchZero or the solving process of CAPD will fail. 

7.2 Bouncing hall 

Example 5.3 can be verified using the implementation by limiting the initial value 
of X 2 to a small interval of width at most 0.01. We verified the model with three 
configurations: 1) with setting a point initial value to X 2 ] 2) with an interval initial 
value of width 0.01; and 3) with a point initial value and the property in which 
the time bound for the G operator was set to [0,100]. Each experiment was run 
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(a) Considered model 


(b) Our method 


(c) dReal 


Figure 4. Results of boundary detection of the bouncing ball example 

1000 times with the initial values of X 2 randomly picked within [2, 7]. The results 
are shown in Table 1. In the table, the second column represents the width of the 

Table 1 

Experimental results (bouncing ball) 


‘f 

width 

# valid 

# unsat 

# unknown 

^ errors 

time 

G[o,io]F[o,5]2 —<0 

0 

330 

670 

0 

0 

0 .1s 

G[o,io]F[o,5]2-xi<0 

0.01 

87 

10 

903 

903 

0 .1s 

G[o,ioo]F[o,5]2 —<0 

0 

186 

20 

794 

794 

0.5s 


initial values. For each r G {valid, unsat, unknown}, the column r” represents 
the number of runs that resulted in r. The column “ 7 ^ errors” represents that how 
many of unknown results are caused by an error in the SearchZero procedure. The 
column “time” shows average timings. 

The computation of each experiment was quite efficient. 

Regarding the rate of inconclusive runs in each experiment, all the verifications 
succeeded in the first experiment. Because we set the point initial values and the 
bounded simulation time, considered trajectories were always enclosed with tight 
intervals, and thus the verification process succeeded even in a situation that was 
close to singular. The result also implied that we did not meet a zeno behavior. 
Contrastingly, around 90% and 80% of the runs in the second and third experiments, 
respectively, resulted in errors. Our verification process with the interval Newton 
failed more often if a trajectory and a guard approached or they became close to 
tangent. Because the model was chaotic, a coarser enclosure of states or a longer 
simulation time increased the possibility to meet such situations. 

In the second and third experiments, most of the successful runs resulted in valid. 
These runs became stable (i.e., there are less difference between the continuous 
trajectories in each step) in the later steps, and thus satisfied the property ip. We 
confirmed that most of the inconclusive runs fell into zeno behaviors. 

It was quite rare that the result of Analintervals became unknown because there 
were always a few (or no) tight boundaries in and they rarely contained zero. 

Next, we experimented with dReal (version 2.14.08), a solver for d-weakened 
SMT problems, for comparison. We consider a portion of the model of (another 
instance of) the bouncing ball such that the initial state is ( 1 . 1 , 0 , 0 ) and the tra¬ 
jectories of the ball and the table become close to tangent (Figure 4 (a)). Next, 
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Figure 5. A trajectory of ATM (left) and the BLTL property to verify (right) 


we analyzed this model by simulating the underlying HA with our method and by 
solving the SMT problem with dReal, respectively. Figure 4 (b) and (c) show the 
computed witness trajectories. Our implementation verified the occurrence of the 
first contact with the floor. dReal computed two enclosures for the possible trajec¬ 
tories of the model; they seemed corresponding to the first and last intersections 
with the guard. The second witness seemed wrong because the guard condition 
became J-sat around t' = 0.5 with 5 = 0.001. The computation of a boundary can 
be troublesome in this way, without the unique existence verification process. 


7.3 Air Traffic Maneuver 

We performed a verification of a simplified model of an air traffic maneuver 
(ATM) [22] in which the number of aircrafts was parameterized. A trajectory of 
(x*,?/*) when m = 4 is illustrated in Figure 5. We verified the following property 
shown in Figure 5. The first line describes that the distance between each pair of 
aircrafts is larger than the threshold 8/m during ||(/?|| = 30 time units. The follow¬ 
ing of the property describes that all aircrafts reach within the circle with radius 5 
within the time interval [0,10], stay there at least 10 time units, and reach outside 
the circle with radius 10 after another 10 time units. We verified 10 times for each 
instance with m = 2, 4, 6, 8. In each verification, we randomly picked a point initial 
value. All runs resulted in valid. The specification of the instances and the results 
are shown in Table 2. The columns vars” and APs” represent the numbers 


Table 2 

Experimental results (ATM) 


m 

# vars 

# APs 

time 

m 

vars 

# APs 

time 

2 

10 

5 

0.5s 

6 

26 

27 

28s 

4 

18 

14 

5.5s 

8 

34 

44 

98s 


of variables in HA and APs in the property. The average timings rose exponentially 
as m increased. 
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8 Conclusions 

We have presented a sound BLTL validation method that assures that all initialized 
trajectories satisfy the property. The proposed method is able to detect a witness 
trajectory that is verified its unique existence with an interval-based ODE inte¬ 
gration and an interval Newton method. We consider the experimental results are 
promising for the practical use. 

In future work, our method and implementation should be improved to allow 
large and uncertain initial values. Examples in a realistic setting should be demon¬ 
strated with the implementation. 
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